home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / newmat.exe / NEWMAT4.CXX < prev    next >
C/C++ Source or Header  |  1991-07-30  |  11KB  |  347 lines
  1. //$$ newmat4.cxx       Constructors, ReDimension, basic utilities
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #include "include.hxx"
  6.  
  7. #include "boolean.hxx"
  8.  
  9. typedef double real;                 // all floating point double
  10.  
  11. #include "newmat.hxx"
  12.  
  13. #include "newmatrc.hxx"
  14.  
  15. //#define REPORT { static ExeCounter ExeCount(__LINE__,4); ExeCount++; }
  16.  
  17. #define REPORT {}
  18.  
  19. //#define REPORT1 { static ExeCounter ExeCount(__LINE__,4); ExeCount++; }
  20.  
  21. // REPORT1 constructors only - doesn't work in turbo and Borland C++
  22.  
  23. #define REPORT1 {}
  24.  
  25. //#define MONITOR(what,storage,store) \
  26. //   { cout << what << " " << storage << " at " << (long)store << "\n"; }
  27.  
  28. #define MONITOR(what,store,storage) {}
  29.  
  30. /*************************** general utilities *************************/
  31.  
  32. static int tristore(int n)                      // els in triangular matrix
  33. { return (n*(n+1))/2; }
  34.  
  35.  
  36. /****************************** constructors ***************************/
  37.  
  38. GeneralMatrix::GeneralMatrix()
  39. { store=0; storage=0; nrows=0; ncols=0; tag=-1; }
  40.  
  41. GeneralMatrix::GeneralMatrix(int s)
  42. {
  43.    REPORT1
  44.    storage=s; tag=-1;
  45.    store = new real [storage]; MatrixErrorNoSpace(store);
  46.    MONITOR("Make (GenMatrix)",storage,store)
  47. }
  48.  
  49. Matrix::Matrix(int m, int n) : GeneralMatrix(m*n)
  50. { REPORT1 nrows=m; ncols=n; }
  51.  
  52. SymmetricMatrix::SymmetricMatrix(int n) : GeneralMatrix(tristore(n))
  53. { REPORT1 nrows=n; ncols=n; }
  54.  
  55. UpperTriangularMatrix::UpperTriangularMatrix(int n)
  56.    : GeneralMatrix(tristore(n))
  57. { REPORT1 nrows=n; ncols=n; }
  58.  
  59. LowerTriangularMatrix::LowerTriangularMatrix(int n)
  60.    : GeneralMatrix(tristore(n))
  61. { REPORT1 nrows=n; ncols=n; }
  62.  
  63. DiagonalMatrix::DiagonalMatrix(int m) : GeneralMatrix(m)
  64. { REPORT1 nrows=m; ncols=m; }
  65.  
  66. Matrix::Matrix(BaseMatrix& M)
  67. {
  68.    REPORT1 CheckConversion(M);
  69.    GeneralMatrix* gmx=M.Evaluate(MatrixType::Rect); GetMatrix(gmx);
  70. }
  71.  
  72. RowVector::RowVector(BaseMatrix& M) : Matrix(M)
  73. {
  74.    if (nrows!=1) MatrixError("Illegal conversion to row vector"); }
  75.  
  76. ColumnVector::ColumnVector(BaseMatrix& M) : Matrix(M)
  77. {
  78.    if (ncols!=1) MatrixError("Illegal conversion to column vector"); }
  79.  
  80. SymmetricMatrix::SymmetricMatrix(BaseMatrix& M)
  81. {
  82.    REPORT1 CheckConversion(M);
  83.    GeneralMatrix* gmx=M.Evaluate(MatrixType::Sym); GetMatrix(gmx);
  84. }
  85.  
  86. UpperTriangularMatrix::UpperTriangularMatrix(BaseMatrix& M)
  87. {
  88.    REPORT1 CheckConversion(M);
  89.    GeneralMatrix* gmx=M.Evaluate(MatrixType::UT); GetMatrix(gmx);
  90. }
  91.  
  92. LowerTriangularMatrix::LowerTriangularMatrix(BaseMatrix& M)
  93. {
  94.    REPORT1 CheckConversion(M);
  95.    GeneralMatrix* gmx=M.Evaluate(MatrixType::LT); GetMatrix(gmx);
  96. }
  97.  
  98. DiagonalMatrix::DiagonalMatrix(BaseMatrix& M)
  99. {
  100.    REPORT1 CheckConversion(M);
  101.    GeneralMatrix* gmx=M.Evaluate(MatrixType::Diag); GetMatrix(gmx);
  102. }
  103.  
  104. GeneralMatrix::~GeneralMatrix()
  105. {
  106.    if (store)
  107.    {
  108.       MONITOR("Free (GenMatrix)",storage,store)
  109.       delete [storage] store;
  110.    }
  111. }
  112.  
  113. CroutMatrix::CroutMatrix(BaseMatrix& m)
  114. {
  115.    REPORT1
  116.    GeneralMatrix* gm = m.Evaluate(MatrixType::Rect); GetMatrix(gm);
  117.    if (nrows!=ncols) MatrixError("Matrix not square");
  118.    d=TRUE; indx=new int [nrows]; MatrixErrorNoSpace(indx);
  119.    ludcmp();
  120. }
  121.  
  122. /**************************** ReDimension matrices ***************************/
  123.  
  124. void GeneralMatrix::ReDimension(int nr, int nc, int s)
  125. {
  126.    REPORT 
  127.    if (store)
  128.    {
  129.       MONITOR("Free (ReDimensi)",storage,store)
  130.       delete [storage] store;
  131.    }
  132.    storage=s; nrows=nr; ncols=nc;
  133.    store = new real [storage]; MatrixErrorNoSpace(store);
  134.    MONITOR("Make (ReDimensi)",storage,store)
  135. }
  136.  
  137. void Matrix::ReDimension(int nr, int nc)
  138. { REPORT GeneralMatrix::ReDimension(nr,nc,nr*nc); }
  139.  
  140. void SymmetricMatrix::ReDimension(int nr)
  141. { REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
  142.  
  143. void UpperTriangularMatrix::ReDimension(int nr)
  144. { REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
  145.  
  146. void LowerTriangularMatrix::ReDimension(int nr)
  147. { REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
  148.  
  149. void DiagonalMatrix::ReDimension(int nr)
  150. { REPORT GeneralMatrix::ReDimension(nr,nr,nr); }
  151.  
  152. void RowVector::ReDimension(int nc)
  153. { REPORT GeneralMatrix::ReDimension(1,nc,nc); }
  154.  
  155. void ColumnVector::ReDimension(int nr)
  156. { REPORT GeneralMatrix::ReDimension(nr,1,nr); }
  157.  
  158.  
  159.  
  160. /********************* manipulate types, storage **************************/
  161.  
  162. int GeneralMatrix::search(const GeneralMatrix* s) const
  163. { REPORT return (s==this) ? 1 : 0; }
  164.  
  165. int AddedMatrix::search(const GeneralMatrix* s) const
  166. { REPORT return bm1->search(s) + bm2->search(s); }
  167.  
  168. int ShiftedMatrix::search(const GeneralMatrix* s) const
  169. { REPORT return bm->search(s); }
  170.  
  171. int NegatedMatrix::search(const GeneralMatrix* s) const
  172. { REPORT return bm->search(s); }
  173.  
  174. int ConstMatrix::search(const GeneralMatrix* s) const
  175. { REPORT return (s==cgm) ? 1 : 0; }
  176.  
  177.  
  178. MatrixType AddedMatrix::Type() const
  179. { REPORT return bm1->Type() + bm2->Type(); }
  180.  
  181. MatrixType MultipliedMatrix::Type() const
  182. { REPORT return bm1->Type() * bm2->Type(); }
  183.  
  184. MatrixType SolvedMatrix::Type() const
  185. { REPORT return bm1->Type().i() * bm2->Type(); }
  186.  
  187. MatrixType SubtractedMatrix::Type() const
  188. { REPORT return bm1->Type() - bm2->Type(); }
  189.  
  190. MatrixType ShiftedMatrix::Type() const
  191. { REPORT MatrixType mteqel(MatrixType::EqEl); return bm->Type() + mteqel; }
  192.  
  193. MatrixType ScaledMatrix::Type() const { REPORT return bm->Type(); }
  194. MatrixType TransposedMatrix::Type() const { REPORT return bm->Type().t(); }
  195. MatrixType NegatedMatrix::Type() const { REPORT return bm->Type(); }
  196. MatrixType InvertedMatrix::Type() const { REPORT return bm->Type().i(); }
  197. MatrixType ConstMatrix::Type() const { REPORT return cgm->Type(); }
  198.  
  199.  
  200. int AddedMatrix::NrowsV() const  { return bm1->NrowsV(); }
  201. int ShiftedMatrix::NrowsV() const { return bm->NrowsV(); }
  202. int TransposedMatrix::NrowsV() const { return bm->NcolsV(); }
  203. int NegatedMatrix::NrowsV() const { return bm->NrowsV(); }
  204. int ColedMatrix::NrowsV() const { return bm->NrowsV() * bm->NcolsV(); }
  205. int DiagedMatrix::NrowsV() const { return bm->NrowsV() * bm->NcolsV(); }
  206. int ConstMatrix::NrowsV() const { return cgm->Nrows(); }
  207.  
  208. int AddedMatrix::NcolsV() const  { return bm2->NcolsV(); }
  209. int ShiftedMatrix::NcolsV() const { return bm->NcolsV(); }
  210. int TransposedMatrix::NcolsV() const { return bm->NrowsV(); }
  211. int NegatedMatrix::NcolsV() const { return bm->NcolsV(); }
  212. int RowedMatrix::NcolsV() const { return bm->NrowsV() * bm->NcolsV(); }
  213. int DiagedMatrix::NcolsV() const { return bm->NrowsV() * bm->NcolsV(); }
  214. int ConstMatrix::NcolsV() const { return cgm->Ncols(); }
  215.  
  216.  
  217. //  Rules regarding tDelete, reuse, GetStore
  218. //    All matrices processed during expression evaluation must be subject
  219. //    to exactly one of reuse(), tDelete(), GetStore() or BorrowStore().
  220. //    If reuse returns TRUE the matrix must be reused.
  221. //    GetMatrix(gm) always calls gm->GetStore()
  222. //    gm->Evaluate obeys rules
  223. //    bm->Evaluate obeys rules for matrices in bm structure
  224.  
  225. void GeneralMatrix::tDelete()
  226. {
  227.    if (tag<0)
  228.    {
  229.       if (tag<-1) { REPORT store=0; delete this; return; }  // borrowed
  230.       else { REPORT return; }
  231.    }
  232.    if (tag==1)
  233.    {
  234.       REPORT  MONITOR("Free   (tDelete)",storage,store)
  235.       if (store) delete [storage] store; store=0; tag=-1; return;
  236.    }
  237.    if (tag==0) { REPORT delete this; return; }
  238.    REPORT tag--; return;
  239. }
  240.  
  241. BOOL GeneralMatrix::reuse()
  242. {
  243.    if (tag<-1)
  244.    {
  245.       REPORT
  246.       real* s = new real [storage]; MatrixErrorNoSpace(s);
  247.       MONITOR("Make     (reuse)",storage,s)
  248.       real* s1=store; real* s2=s; int i=storage;
  249.       while(i--) *s2++ = *s1++;
  250.       store=s; tag=0; return TRUE;
  251.    }
  252.    if (tag<0) { REPORT return FALSE; }
  253.    if (tag<=1)  { REPORT return TRUE; }
  254.    REPORT tag--; return FALSE;
  255. }
  256.  
  257. real* GeneralMatrix::GetStore()
  258. {
  259.    if (tag<0 || tag>1)
  260.    {
  261.       real* s = new real [storage]; MatrixErrorNoSpace(s);
  262.       MONITOR("Make  (GetStore)",storage,s)
  263.       real* s1=store; real* s2=s; int i=storage;
  264.       while(i--) *s2++ = *s1++;
  265.       if (tag>1) { REPORT tag--; }
  266.       else if (tag < -1) { REPORT store=0; delete this; } // b